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ABSTRACT  (Cow Him*  m  nvarn  MM  if  aumwr  m4  IMawlliy  Bp  Bleak  mwAer) 

_  Magnetohydrodynamic  stability  of  force  free  spheromak  plasmas  to  n*l  free  boundary 
tilt  and  radial  drift  modes  is  studied.  The  exterior  region  of  open  field  lines  is  treated 
either  aa  a  conducting  plasma  or  as  a  vacuum.  In  the  former  case  the  effect  of  line  tying 
of  the  open  field  lines  is  present  and  is  found  to  improve  stability  greatly.  The  equilibria 
which  have  optimum  stability  properties  to  the  tilt  and  shift  modes  have  s  length  to 
dtaaMtsr  ratio  of  0.6,  and  can  be  stable  with  conducting  walla  at  a  relatively  large  dittriir^t. 
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TILT  AND  SHIFT  MODE  STABILITY  WITH  LINE  TYING 


I .  INTRODUCTION 

A  spheromak  is  a  compact  torus  magnetic  configuration 
whose  toroidal  and  poloidal  magnetic  fields  are  produced 
by  plasma  currents  confined  by  an  external  vertical  field. ^ 
The  advantages  it  shares  with  other  compact  torus  designs 
are  compactness  and  simplicity  of  design  of  external  coils, 
walls  and  plasma,  as  well  as  a  natural  divertor.  A  plot 
of  the  flux  surfaces  of  a  typical  spheromak  configuration 
is  shown  in  Fig.  1.  Note  the  two  classes  of  field  lines, 
those  on  closed  flux  surfaces  and  the  open  field  lines, 
separated  by  a  surface  we  call  the  separatrix. 

It  has  been  found  that  a  spheromak  can  be  stable  to 
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internal  modes  if  it  is  oblate.  Even  an  oblate  spheromak 
requires  an  external  conductor  to  stabilize  the  free¬ 
boundary  modes.  The  question  of  how  close  the  conducting 
shell  must  be  and  what  shape  it  is  allowed  to  take  has  an 
important  bearing  on  spheromak  experiments,  and  on  the 
viability  of  proposed  spheromak  reactor  designs.  The 
moving  ring  field  reversed  reactor,*  for  example,  would 
require  a  cylindrical  shell,  with  the  ends  of  the  cylinder 
either  far  away  or  tfbsent  altogether.  More  conventional 
reactor  designs  would  allow  a  shell  which  almost  completely 
surrounds  the  plasma,  but  would  nonetheless  place  constraints 
on  how  far  the  shell  must  be  from  the  fusion  plasma. 
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In  this  paper  we  numerically  study  wall  stabilization 

of  the  external  n  ■  1  nagnetohydrodynamic  (MHD)  inodes  of  the  spheromak, 

with  the  effect  of  line  tying  due  to  plasma  on  the  open  field 

lines.  Previous  work  indicates  that  the  n  =  1  modes  are 

the  most  difficult  of  the  external  ideal  MHD  inodes  to 

stabilize.*  The  n  ■  1  tilt^'®'^  and  radial  shift^  modes 

have  also  been  seen  in  spheromak  experiments. 

The  first  studies  of  MHD  stability  of  spheromaks, 

using  a  nearly  spherical  equilibrium  with  7  x  ■  uJJ  and 

2 

u  constant,  are  due  to  Rosenbluth  and  Bussac.  They 

found  that  slightly  prolate  equilibria  are  unstable  to  an 

internal  tilting  mode,  while  slightly  oblate  equilibria 

require  a  shell  of  radius  R(l-0.2e)  to  stabilize  the 

external  tilting  mode,  where  R  is  the  plasma  radius  and 

e  is  the  ellipticity  of  the  elliptical  plasma  surface. 
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Later  work  by  Finn  et  al.  and  by  Bondeson  et  al.  has 

shown  that  in  a  cylindrical  conducting  wall  with  endplates 
the  internal  tilt  mode  is  unstable  as  long  as  the  ratio 
of  length  to  radius  exceeds  1.67. 

The  existence  of  a  q  ■  1  surface  in  the  plasma  (q 
is  the  usual  safety  factor)  has  been  shown  to  be  a  destabi¬ 
lizing  factor  in  oblate  plasmas.*  The  equilibria  studied 
in  the  present  paper  have  q  less  than  unity  at  the  O-point 
and  q  decreasing  toward  the  plasma  edge,  so  that  no  such 
q  ■  1  surface  exists.  Our  equilibria  are  force-free. 
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Tilting  modes  have  been  observed  in  several  experi¬ 
ments  to  date,  notably  on  the  PS-1  device  at  the  University 
of  Maryland^  and  on  plasmas  produced  by  magnetized  coaxial 
plasma  guns  at  Los  Alamos**  and  Livermore.^  On  the  PS-1 
experiment  another  n  =  1  instability  has  also  been  observed, 
a  radial  shift  mode.^  (This  should  be  distinguished  from 
the  axial  shift,  which  is  axi symmetric . ) 

The  study  described  in  this  paper  extends  that  of 
Rosenbluth  and  Bussac  in  considering  more  general  separatrix 
shapes  and  in  including  two  additional  physical  effects 
that  turn  out  to  be  quite  stabilizing.  The  first  effect 
is  that  we  consider  equilibria  whose  current  goes  to  zero 
smoothly  at  the  separatrix,  or  even  at  an  inner  flux 
surface  (so  that  there  is  a  "flux  hole");  i.e.,  u  *  p(*). 

This  smoothing  of  the  current  profile,  which  is  almost 
certainly  present  in  experiments,  has  an  important  stabi¬ 
lizing  effect.  This  is  plausible  because  a  large  part  of 
the  free  energy  for  the  free  boundary  tilt  is  due  to  the 
magnetic  pressure  imbalance  at  the  perturbed  X-point. 

The  second  physical  effect  we  include  is  line  tying 
of  the  open  field  lines.  That  is,  we  consider  the  external 
region  to  be  filled  with  conducting  plasma  (carrying  zero 
equilibrium  current)  rather  than  vacuum.  Since  there  is 
a  normal  component  of  magnetic  field  at  the  conducting 
walls  (see  Fig.  1),  this  effect  is  stabilizing.  Note  that 
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this  configuration  avoids  a  serious  problem  of  conventional 
line-tied  systems,  that  impurities  can  flow  into  the  hot 
plasma  from  the  metal  wall.  Although  the  closed  field 
lines  on  which  the  equilibrium  current  flows  are  not  them¬ 
selves  line-tied,  we  find  that  line  tying  of  the  open  field 
lines  greatly  improves  stability  to  n  =  1  modes. 

Recently,  other  groups  have  also  performed  stability 
computations  on  the  tilt  and  radial  shift  modes.10,11  The 
equilibrium  models  permitted  a  flux  hole  but  did  not  permit 
finite  current  arbitrarily  close  to  the  separatrix.  But 
the  most  important  difference  between  this  work  and  the 
work  we  present  in  this  paper  is  our  inclusion  of  the 
stabilizing  effect  of  line  tying  of  the  open  field  lines. 

Our  results  without  line  tying  (i.e.,  with  a  vacuum 
exterior)  are  consistent  with  the  results  of  Refs.  10  and  11. 
(Exact  comparison  is  difficult  due  to  differences  in  pro¬ 
files  and  wall  shapes.)  However,  we  find  that  the  equi¬ 
libria  are  considerably  more  stable  with  line  tying.  For 
example,  one  of  our  equilibria  optimized  for  tilt  and  shift 
modes  is  marginally  stable  without  line  tying  when  the  wall 
radius  is  1.4  times  the  separatrix  radius  and  the  spacing 
between  the  axial  walls  is  1.1  times  the  separatrix  length. 
With  line  tying,  the  last  figure  is  1.5. 

In  Section  II  we  describe  the  equilibrium  model, 
in  particular  the  profiles  and  boundary  conditions  used. 
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In  Section  III  we  describe  modified  linear  magneto¬ 
hydrodynamic  equations  for  force  free  plasmas  and  the  time 
dependent  code  FFMHD  that  we  use  to  integrate  these  equa¬ 
tions  of  motion.  We  describe  the  special  properties  of 
the  modified  magnetohydrodynamic  equations  that  make  them 
particularly  suitable  to  a  system,  such  as  the  spheromak, 
that  has  magnetic  neutral  points  (the  X  points  on  the 
symmetry  axis) .  We  also  describe  how  the  exterior  region 
may  be  conveniently  treated  as  a  perfectly  conducting 
plasma  or  a  vacuum  in  this  framework. 

In  Section  IV  we  discuss  the  physical  properties  of 
the  tilt  and  radial  shift  modes;  in  particular  we  show 
results  from  FFMHD  pertaining  to  the  stabilizing  effect 
on  the  tilt  mode  of  smoothing  the  current  near  the 
separatrix. 

In  Section  V  we  describe  in  detail  results  obtained 
using  FFMHD  on  the  effect  of  elongation  of  the  separatrix 
on  the  tilt  and  shift  modes.  There  is  an  optimum  elonga¬ 
tion  (separatrix  half  length  to  radius  ratio),  equal  to 
about  0.6,  nearly  independent  of  other  parameters.  The 
optimum  elongation  without  line  tying  is  approximately 
the  same,  but  the  walls  must  be  much  closer  to  achieve 
stability.  We  also  present  marginal  stability  results  in 
which  we  vary  the  elongation  of  the  cylindrical  wall.  We 
find  that  line  tying  allows  the  axial  wall  to  be  removed 
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to  a  large  distance,  as  required  by  the  moving  ring  reactor 
scheme.  We  conclude  that  line  tying  in  the  radial  wall  is 
particularly  stabilizing. 


I 
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II.  EQUILIBRIUM  MODEL 


Since  the  instabilities  under  consideration  are  current 
driven,  we  restrict  our  attention  to  force  free  equilibria, 
i.e. ,  those  equilibria  having  j  =  yB.  However,  we  generalize 
the  equilibria  considered  in  Refs.  2,  8,  and  9  by  not  requiring 
u  to  be  a  constant  in  the  plasma  and  by  not  requiring  the 
separatrix  to  be  nearly  spherical.  Instead,  we  keep  only  the 
obvious  requirements  that  u  be  constant  on  flux  surfaces  and 
go  to  zero  continuously  as  the  separatrix  is  approached  from 
the  inside. 

From  the  general  axisymmetric  representation  in  cylindrical 
coordinates  (r,  0,  z) 

B  =  7’4j  x  V6  +  g(t)V0  (1) 

% 

we  conclude 


j  =  V  x  B 


-  -  A*i|>?0  +  g'(*)7ii> 

x  78  , 

(2) 

where  primes 

denote  differentiation 

with  respect  to 

'l>  =  rAg,  the 

poloidal  flux,  and  A*y 

=  r27.  (r~27’p)  . 

The 

requirement  ' 

j  =  'u('ii)  B  gives 

u  (  v)  -  g'  (v)  , 

(3a) 

=  -  g(y)  g'(p) 

. 

(3b) 

8 


( 
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The  last  equation  is  a  specialization  of  the  Grad-Shaf ranov 
equation  to  force  free  fields.  We  choose  j  where 
i^<0  is  the  magnitude  of  the  flux  hole. 


for  J'O,  (4a) 


=  0 


for  J>0,  (4b) 


where  '^q<0  is  the  minimum  of  ^ ,  which  occurs  at  the 
magnetic  axis  or  0-point;  this  model  gives 


u  (*) 


aBeo  _ $ _ 

vo  [u/*0)2  +  s2]** 


for  ;<0,  (5a) 


=  0  for  i>0.  (5b) 

The  boundary  conditions  are  V;/rZ  =  0  at  z  =  0  (reflection 
symmetry),  J  =  0  at  r  =  0  (from  ;■  =  r  A. )  and  the  value  of 
'l>  as  a  function  of  z  at  r  =  a  (the  radial  wall),  and  as  a 

function  of  r  at  z  =  L  (the  endplate) .  The  separatrix  is 

given  by  <i>  -  0.  From  (5)  we  see  that  the  equilibrium 
current  is  zero  on  the  open  field  lines  (v  >  0) ,  that  u 
is  nearly  constant  near  ^  =  Pq  and  that  u  is  linear  in  .■ 

near  the  flux  hole  [(^h  -  <M/|'vq|  £  5]*  Notice  that  if 

5  =  0,  =  0,  the  poloidal  current  is  discontinuous  at 

the  separatrix.  The  specification  of  g  as  a  function  of 
the  ratio  P/pq  as  in  (4)  provides  convergence  of  the  itera¬ 
tive  scheme  used  to  solve  (3b).  See  Ref.  12. 
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One  set  of  boundary  conditions  at  r  *  a  and  z  ■  L, 
expressed  in  terms  of  dimensionless  variables  a  *  ;/(*jBga2) 
where  Bg  is  a  reference  field,  is 


'P 


,,  ,  .  -3  .2  2. -3/2 

(l-.6e)r0  -  (r  +  z  ) 


+  e(4z2  -  r2)[-.175r'5  -  .425r2(r‘ 


(6) 


2 

z 


I 

r 


For  |e|<<l,  these  boundary  conditions,  for  :  *  0,  give  the 

nearly  spherical  constant  u  equilibria  of  Rosenbluth  and 

2  2 
Bussac,  with  plasma  surface  given  by  r  =  r^tl  +  cos  ,), 

where  x  =  tan  ^(r/z).  Equilibria  with  •  >  0  are  prolate, 

and  those  with  e  <  0  are  oblate.  Another  set  of  boundary 

conditions  that  we  use  is 


=  CLr2(l 


+  C2z2) 


(7) 


For  this  set  of  boundary  conditions,  equilibria  with  small 

have  a  separatrix  which  is  near  the  wall.  In  fact,  for 

=  i  =  0 ,  the  equilibria  are  those  obtained  analytically 
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by  Finn,  Manheimer  and  Ott  and  by  Pondeson  et  al.  For 
>  0,  the  parameter  C 2  controls  the  oblateness  of  the 
equilibrium  by  specifying  the  amount  of  flux  passing 
through  the  radial  wall  r  =  a  between  z  3  0  and  z  ■  L. 


The  flux  surfaces  u  =  const,  of  a  typical  equilibrium 
are  shown  in  Fig.  1.  Notice  that  the  separatrix  (separating 


open  field  lines  and  those  on  closed  flux  surfaces)  intersects 
the  z-axis  at  magnetic  neutral  points,  which  we  call  X-p  'ints, 
one  on  each  side  of  the  plane  of  sysasetry  z  -  0. 

Our  model  (4)  produces  plasmas  with  fairly  flat 
toroidal  current  profile.  We  observe  that  for  all  plasmas 
oblate  enough  to  be  stable  to  the  internal  tilting  mode, 
the  safety  factor  q  is  l<*  i  than  unity  throughout  the  plasma. 
Therefore,  for  n  -  1  modes  no  mode  rational  surfaces  exist 
in  the  plasma.  Furthermore,  we  do  not  expect  that  adding 
pressure  with  a  fairly  flat  profile  to  a  plasma  with  fixed 
separatrix  will  cause  the  safety  factor  q  to  }o  above  unity. 
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III.  MODIFIED  MAGNETOHYDRODYNAMIC  EQUATIONS 

In  this  section  we  describe  the  modified  magneto- 
hydrodynamic  equations  of  motion  which  are  integrated  in 
time  by  the  numerical  code  FFMHD  until  the  most  unstable 
mode  dominates.  For  more  detail  on  the  code,  see  Appendix 
A.  The  equation  of  motion  for  a  force  free  equilibrium 
and  ideal  Ohm's  law  are 

_  £  ■  :  j  *  B  ♦  j  ■  \B ,  (8a) 

e  ♦  l  •  B  -  0,  (8b) 

where  :f  is  the  Eulenan  plasma  displacement,  .  is  the 
density  and  dots  denote  derivation  with  respect  to  time. 
From  (8)  and  Max-*  11'  s  equations,  we  conclude  that  the 
perturbed  fields  can  be  represented  by  a  perturbed  vector 
potential  A  ■  £  *  B  and  zero  scalar  potential.  From  (8a), 
therefore,  we  find 

^  A  -  (  ,)  B  -  \]JA,  (9) 

where  i  represents  components  perpendicular  to  the  unper¬ 
turbed  field  B.  This  is  the  basic  equation  of  motion  for 
the  plasma  which  is  integrated  in  time  by  the  numerical 
code  FFMHD.  The  Taylor  stability13  ot  equilibria  with 
uniform  ■„  was  computed  in  Ref.  8  by  integrating  (9)  without 


j 
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the  operation  denoted  by  x;  it  was  shown  that  Taylor 
stability  of  such  equilibria  is  equivalent  to  ideal  mag¬ 
netohydrodynamic  stability  if  no  mode  rational  surfaces 
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exist  in  the  plasma. 

Because  of  the  toroidal  symmetry  of  the  equilibrium 
(V>d  *  0),  the  cylindrical  components  of  £  or  of  iA  may 
be  assumed  to  behave  as  exp(inS),  where  n  is  called  the 
toroidal  mode  number. 

For  external  or  free  boundary  modes,  integrating 
(9)  to  obtain  :A  rather  than  (8a)  to  obtain  the  displace- 

ment  \  has  a  further  advantage  over  those  mentioned  in 

■> 

Ref.  8  (for  internal  modes).  As  discussed  in  Ref.  2,  a 
slightly  oblate  equilibrium  can  be  unstable  to  tilting 
because  the  internal  tilt  motion  is  nearly  marginally 
stable,  and  this  motion  produces  a  magnetic  pressure 
imba’ance  across  the  perturbed  X-points.  This  additional 
source  of  free  energy  drives  the  instability.  Therefore, 
we  might  expect  fairly  singular  behavior  of  ;  near  the 
X-points,  since  5  should  be  large  in  that  vicinity,  but 
must  be  zero  on  the  z-axis.  (The  latter  property,  which 
must  hold  for  every  vector  with  toroidal  mode  number 
n  ^  1,  is  discussed  further  in  Appendix  A.)  The  displace¬ 
ment  plots  shown  in  Fig.  2  show  that  this  singular  behavior 
is  indeed  found,  and  therefore  that  numerical  problems 
would  almost  certainly  be  encountered  in  attempting  to 
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integrate  (da) .  On  the  other  hand,  6A  is  well  behaved 
near  the  X-points,  and  no  such  numerical  problems  are 
observed  in  integrating  (9) .  This  property  and  the  fact 
that  no  mode  rational  surfaces  exist  in  the  plasma  allows 
us  to  use  a  simple  cylindrical  (r,z)  grid  rather  than  one 
matched  to  the  flux  surfaces  of  the  equilibrium. 

We  have  two  optional  ways  of  treating  the  exterior 
region  of  open  field  lines.  The  first,  called  the  line 
tying  option,  is  to  treat  the  region  as  a  perfectly  con¬ 
ducting  plasma  carrying  no  equilibrium  current.  That  is, 
since  u(y)  is  zero  on  the  open  field  lines  [c.f.  equation 
(5,] ,  Wl  merely  integrate  (9)  over  the  whole  region  without 
regard  to  the  location  of  the  separatrix.  The  boundary 
conditions  ft  *  SE  *  0,  n  •  ;  ■  0  at  the  conducting  impene- 
trable  wall,  together  with  the  fact  that  n  •  B  is  not  zero 
imply  £  *  0  at  the  wall;  this  with  the  condition  B  •  cA  *  0 
yields  the  boundary  condition  on  the  vector  potential 
JA  *  0.  Using  this  boundary  condition  it  is  easy  to  see 
that  (9)  has  the  usual  energy  principle,  namely  that  the 
normal  modes  are  extrema  of 


where 


5W  -  [cB2  -  u(v»)5A  •  6e]d3x  , 
T  -  jj*(D/B2)5A2  d3x 


(10) 

(11a) 

(lib) 
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d. 


and  where  6A  satisfies  the  constraint  B  •  6A  *  0.  Using 

X  A,  'V, 

the  toroidal  symmetry  of  the  equilibrium,  we  can  write  the 
toroidal  dependence  explicitly,  i.e.,  <5A  *  cos  n8  - 
5A1  sin  n0 ,  -iB  a  5Br  cos  n0  -  5B1  sin  n0.  Using  these  real 

\i  \  \  'V, 

variables  in  (10),  (11)  is  equivalent  to  using  the  complex 
variables  6A  =»  (5Ar  +  ifiA1)  exp(ind) ,  6B  =  (5Br  +  ifiB1) 

<v,  %  'V  %  'V 

exp(in0)  with  (10)  but  with 


5W 


V  (*i») 


Re  ( 6A  * 


d3x 


T 


kfi: 


/B2)  |  5A  | 


d3x 
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(12a) 


and  B  •  5  A  »  0. 

The  exterior  plasma  model  includes  the  important 
stabilizing  effect  of  line-tying  of  the  open  field  lines. 
One  way  to  see  that  this  is  a  stabilizing  influence  is  to 
note  that  having  conducting  plasma  in  the  exterior  region 
requires  that  the  ideal  magnetohydrodynamic  flux  constraint 
!B  *  7  *  (£  x  b)  or  5 E 1 1  *  0  holds  throughout  the  region, 
while  this  is  not  the  case  if  the  exterior  region  is 
vacuum.  Also  recall  that  the  boundary  condition  5At  *  0» 
together  with  (8a)  and  ^  =  p B  implies  that  the  total  plasma 
displacement  E,  is  zero  at  the  wall.  (If  the  wall  is  a 
flux  surface,  only  the  normal  component  of  £  must  be  zero.) 
As  we  shall  discuss  further,  the  ’boundary  conditions  when 
the  conducting  wall  surrounds  a  vacuum  are  A  *  eA  *  0, 
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regardless  of  whether  a  normal  component  of  the  equilibrium 
field  exists.  In  Appendix  B  we  establish  more  rigorously 
the  stabilizing  effect  of  line  tying. 

Various  different  normalizations,  i.e.,  models  for 
the  density  p(r,z)  are  possible.  The  choice  affects  growth 
rates  but  not,  of  course,  marginal  stability.  The  simplest 

such  model  has  p  =  const.  Another  model  which  is  ver; 

2  2 

useful  has  B  /p  ~  VA  “  const,  (the  Alfvdn  speed).  Whereas 
this  model  has  the  unphysical  effect  of  allowing  equilibrium 
density  variations  along  the  field  lines,  it  has  the  advan¬ 
tage  of  producing  growth  rates  of  the  order  of  three  times 
those  of  the  constant  density  model.  This  advantage  is 
due  to  the  fact  that  the  density  goes  to  zero  at  the  X- 
points,  where  a  large  destabilizing  contribution  occurs. 

Any  model  having  p  =  p('p)  and  with  small  density  near  the 
X-points  would  have  a  very  large  Alfven  speed  v  near  the 
separatrix  and  would  therefore  have  a  very  stringent 
Courant  condition.  For  more  details  see  Appendix  A.  Of 
course,  the  unphysical  nature  of  the  second  model  has  no 
effect  on  the  marginal  stability  computations.  We  expect 
that  the  growth  rates  for  a  system  whose  density  decreases 
monotonically  to  a  finite  value  at  the  separatrix  would  be 
bounded  by  the  growth  rates  of  our  two  extreme  models. 
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The  alternative  model  for  the  exterior  region  is  a 

14 

vacuum.  The  usual  energy  principle  treatment  of  the  vacuum 
i3  to  minimize  the  vacuum  contribution  to  5W, 


=  /<* 


x  cA)2  d3x. 


that  is  satisfy  V  *  7  x  5A  «  0,  subject  to  the  boundary 
conditions  ft  •  5  B  =  -  ft  *  5A  at  the  plasma-vacuum  inter- 

<\j  ^  -v 

face  and  ft  *  6A  =  0  at  the  conducting  wall.  Then,  stability 
% 

is  determined  by  minimizing 


Ui  =  6W/T, 


where  oW  consists  of  the  plasma  contribution  (11a)  plus 

1  f  ?3 

5Wv  and  where  T  is  still  given  by  (lib)  or  ^  /  c£  d  x.  Our 
method,  using  the  time  dependent  code  FFMHD,  is  to  inte¬ 
grate  (9)  on  the  closed  flux  surfaces  and,  on  the  open 
field  lines,  to  integrate 


06 A  =  -  B' 


exp(-pu>)  6 


where  p  is  an  adjustable  parameter  and  ||  denotes  components 
parallel  to  B.  Clearly  the  total  differential  equation  (9) , 
(15)  is  continuous  across  the  separatrix  4/  *  0.  Far  from 
the  separatrix,  specifically  when  the  exponential  factor 
in  (15)  is  negligible,  (15)  reduces  to  the  form  of  the 
equation  of  motion  without  the  constraint  B  •  5A  ■  0.  We 


I 


4. 


generally  pick  p  large  so  that  there  is  only  a  very  small 
region  immediately  outside  the  separatrix  where  the  expo¬ 
nential  factor  in  (15)  is  nonnegligible .  There  are  cases, 
in  fact,  in  which  p  can  be  so  large  that  the  exponential 
factor  is  negligible  at  every  grid  point  in  the  exterior 
region.  However,  occasionally  it  is  necessary  to  use  a  more 
moderate  value  of  p  to  avoid  numerical  instabilities. 

Taking  the  exponential  factor  in  (15)  to  be  zero, 

we  easily  see  that  equations  (9),  (15)  have  an  energy 

2 

principle  represented  by  w  *  5W/T  ,  where  iw  consists 
of  the  same  plasma  and  vacuum  terms  as  in  (14)  but  where 


T' 


* y*(o/B2)6A^  d3x 
P 


+  J  (0/ B2)6A2 

V 


d3x. 


(16) 


Here,  the  notation  serves  to  emphasize  the  fact  that 
B  •  iA  *  0  ir  ;ide  the  plasma,  but  that  -no  such  constraint 
exists  in  the  vacuum.  In  Appendix  B  we  show  that  the  condi¬ 
tion  that  ft  *  5A  be  continuous  across  the  separatrix  (equiv- 
alent  to  ft*  £  B  «  -  ft  *  5A)  is  indeed  satisfied  by  (9) , 

(15)  even  when  exp(-pv)  is  taken  to  be  zero.  The  density 
o  in  the  vacuum  is,  of  course,  fictitious;  however  £W  con¬ 
sists  of  the  same  two  terms  as  the  conventional  energy 
principle  and  therefore  predicts  marginal  stability  correctly. 
In  principle,  this  fictitious  density  could  be  made  small 
in  order  to  compute  actual  growth  rates,  but  the  resulting 
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large  Alfv^n  speed  would  produce  a  very  stringent  Courant 
condition. 

There  are  two  other  features  of  the  code  FFMHD  worth 

discussing.  The  first  is  that  we  generally  add  an  addi- 
2 

tional  term  yrt  p  6a  to  the  right  hand  side  of  the  equations 
u 

of  motion  (9)  in  the  case  of  plasma  exterior  and  (9),  (15) 

in  the  case  of  vacuum  exterior.  Th’ *  modification  speeds 

up  convergence,  i.e.,  the  evolution  to  a  pure  normal  mode, 

if  only  one  unstable  mode  exists.  In  this  case  the  optimum 

value  of  Yq  is  l^!/  w^ere  is  fre<3uency  the 

second  normal  mode,  the  mode  with  the  second  largest 
2 

eigenvalue  ~j>  .  If  two  unstable  normal  modes  exist,  or  if 

the  second  mode  has  ui  ^  0,  this  additional  term  does  not 

help.  Indeed,  we  have  observed  that  when  a  q  =  1  surface 

exists  in  the  plasma,  where  q  is  the  usual  safety  factor, 

the  convergence  rate  to  n  =  1  instabilities  : s  not  aided 

by  having  finite  y^,  presumably  because  continuum  eigen- 

2 

functions  with  u  0  exist.  However,  throughout  this  paper 

we  deal  with  short  (oblate)  plasmas  where  q  <  1  throughout 

the  plasma.  In  addition,  in  cases  for  which  no  unstable 

normal  mode  exists,  using  Yq  enables  us  to  compute  the 

2 

normal  mode  with  smallest  j  . 

The  other  noteworthy  feature  of  the  code  is  that  we 
compute  growth  rates,  for  the  purpose  of  extrapolation  to 
marginal  stability,  by  (10)  rather  than  merely  observing  the 


* 


exponential  growth  rate  directly  after  a  pure  mode  is 
obtained.  We  call  the  former  growth  rate  Yn  and  the  latter 
Yt-  That  is,  our  code  is  essentially  an  energy  principle 
code  where  the  optimum  trial  function  is  obtained  by  inte¬ 
grating  the  modified  equations  of  motion.  Our  expectation 
that  the  growth  rates  Yn  computed  by  this  method  will  be 
more  accurate  is  based  on  Rayleigh's  principle,  namely  that 
the  error  in  the  eigenvalue  computed  by  a  variational  prin¬ 
ciple  is  of  the  order  of  the  square  of  the  error  (e.g., 
grid  errors)  in  the  eigenfunction.  On  the  other  hand,  the 
growth  rate  Yfc  observed  directly  in  the  code  should  have 
errors  of  the  same  order  as  the  eigenfunction.  We  have 
performed  convergence  tests  with  increasingly  finer  grids 

in  r,z  and  we  have  always  found  that  7  is  more  accurate. 

n 

For  more  detail  see  Appendix  A. 

Computing  growth  rates  Yn  by  (10)  also  allows  us  to 
run  for  shorter  times,  since  a  first  order  error  in 
the  eigenfunction  due  to  incomplete  convergence  to  the 
most  unstable  normal  mode  leads  only  to  second  order  errors 
in  the  growth  rate.  This  property  also  holds  for  Yfc  if 
the  growth  rate  is  computed  by  observing  the  e-folding 
rate  of  T  (lib) ,  because  of  orthogonality.  This  property 
does  not  hold  for  any  other  quadratic  form.  However,  the 
errors  in  this  method  of  computing  the  growth  rate  become 
amplified  near  marginal  stability  if  Yq  is  finite.  It  is 
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easily  shown  that  the  error  in  relative  to  Yq  becomes 
first  order  in  the  error  in  the  eigenfunction  near  marginal 
stability.  The  error  in  yfc  relative  to  yfc  itself  is,  of 
course,  even  larger. 
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IV.  TILT  AND  SHIFT  MODES 

In  this  section  we  present  numerical  studies  of  the 

stability  of  free  boundary  n  =  1  modes  with  line  tying  on 

the  open  field  lines.  We  restrict  this  study  to  oblate 

equilibria  because  prolate  equilibria  are  unstable  to 

internal  tilt  modes.  For  the  equilibria  we  consider  (4), 

the  safety  factor  q  is  below  unity  throughout  the  plasma, 

and  therefore  the  rather  singular  behavior  near  the  q  =  1 

4  9 

surface  typical  of  internal  kink  modes  '  does  not  occur. 

As  we  have  discussed  in  a  different  context  i .  •  Sec. 

Ill,  a  free  boundary  tilt  motion  may  be  unstable  with  walls 

fairly  close  to  the  plasma  because  the  nearly  marginal  internal 

tilt  motion  produces  a  magnetic  pressure  imbalance  across 

the  perturbed  X-points.  A  modification  of  the  displacement 

^  to  take  advantage  of  this  free  energy  ran  provide  SW<0, 

2 

i.e.,  instability.  In  fact,  the  method  of  Rosenbluth 
2 

and  Bussac  of  computing  stability  of  equilibria  with 
uniform  u  is  based  on  computing  this  magnetic  pressure 
imbalance.  However,  their  method  does  not  generalize  to 
u(i')  in  any  obvious  way. 

The  poloidal  displacement  £  =  z  of  a  tilt 

mode  in  a  typical  equilibrium  is  shown  in  Figure  2.  (The 

displacement  is  computed  from  the  growing  solution  6A  of 

2 

(9)  by  £  =  B  *  i A/B  .)  Notice  that  the  imaginary  part  of 
£,  i.e.,  the  sin  6  Fourier  component  is  indeed  a  nearly 
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rigid  rotation  together  with  an  inward  motion  at  the 

rotated  X-point.  The  real  part,  i.e.,  the  cos  9  component, 

is  a  nearly  rigid  radial  shift. 

Our  general  procedure  for  studying  stability  for  a 

given  equilibrium  is  to  fix  the  position  of  the  radial  wall 

r  =  rw  and  perform  several  stability  runs  with  the  axial 

wall  z  =  z^  successively  closer.  We  then  extrapolate  to 

2  2 

marginal  stability.  The  square  of  the  growth  rate  y  =  -w 

n 

from  (10)  is  usually  quite  linear  near  marginal  stability. 

Typical  growth  rates  found  for  tilting  modes  with 
2  2 

normalization  v^  =  B  /p  =  const,  have  yL/vA  v  l.  For  B  ^ 

14 

2  kg,  n  ^  10  and  L  ^  50  cm,  as  in  the  Los  Alamos  experiment, 
we  find  y  ^  106  sec  1.  ^s  we  have  discussed  earlier,  with 
the  p  =  const,  normalization,  this  is  lower  by  a  factor  of 
order  3. 

Constant  u  equilibria  of  Ref.  2,  i.e.,  with  (4), 

6=0,  4^  =  0,  with  boundary  conditions  (6)  with  e  <  0 
required  a  tight  fitting  shell.  For  example  for  i  =  -0.4, 
with  the  radial  wall  nearly  touching  the  separatrix  at  the 
midplane  rw  =  rg,  marginal  stability  occurs  when  the  axial 
wall  nearly  touches  the  separatrix  at  the  X-point,  i.e., 
z  =  z  •  This  is  in  qualitative  agreement  with  the  results 

W  5 

of  Ref.  2,  namely  that  a  spherical  wall  must  be  at 
2  2  h 

(r  +  z  )  =  rg(l  +  . 2  j  e  j )  for  marginal  stability.  Of 

course,  qualitative  agreement  is  all  we  expect  since  our 

i 
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wall  is  cylindrical  rather  than  spherical  and  since  <:  *  -0.4 
is  not  small  enough  in  absolute  value  for  the  treatment  of 
Ref.  2  to  be  accurate. 

In  Figure  3  we  show  z ,  the  position  of  the  axial  wall 
required  for  stability  as  a  function  of  <5 ,  the  parameter  of 
(4) ,  (5)  that  determines  the  smoothing  of  the  poloidal  current 
(which  is  discontinuous  if  6  =  0)  near  the  separatrix.  The 
radial  wall  was  nearly  touching  the  separatrix,  rv/rs  =  1.01. 
The  stability,  as  measured  by  zw,  improves  dramatically  with 
•5.  Notice  that  the  point  zg  where  the  separatrix  intersects 
the  z-axis  (at  the  X-points)  decreases  with  5,  because  the 
total  current  decreases.  This  stabilization  due  to  smoothing 
the  current  in  the  neighborhood  of  the  separatrix  is  plausible, 
because  the  destabilizing  force,  namely  the  pressure  imbalance 
at  the  X-points,  is  localized  near  the  separatrix. 

In  Figure  4  we  show  the  axial  wall  position  z required 
for  stability,  for  r^r  *1.01  and  '  *  0.15,  as  a  function  of 
the  elongation  parameter  of  (6).  Note  that  there  is  a 
dramatic  improvement  with  oblateness  that  begins  around 
*:  *  -.15.  Improvement  for  large  is  reasonable,  since 

the  results  of  Ref.  8  show  that  the  energy  of  the  n  =  1 
internally  tilted  equilibrium  continues  to  rise  relative 
to  the  axi symmetric  state  as  the  plasma  length  decreases, 
unt'l  extremely  short  plasmas  are  obtained.  It  is  worth 
mentioning  here  that  there  are  no  instances  in  which  z 
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must  be  less  than  z  ,  i.e.,  the  axial  wall  intersecting  the 

s 

separatrix,  to  stabilize  external  tilting  modes,  even  though 
having  z  slightly  less  than  z  would  still  leave  quite  a 

W  S 

large  exterior  region.  .his  is  probably  due  to  the  fact 
that  the  displacement  of  the  boundary  is  large  in  a  neighbor¬ 
hood  of  the  X-point  on  the  z-e  is  This  feature  is  illus¬ 
trated  in  Figure  4,  where  z w  =  Zg  is  required  for  marginal 
stability  over  a  fairly  large  range  of  e,  -0.'5  <  e  <  0, 
even  though  zg  is  decreasing  with  |  e  |  . 

The  results  shown  in  Figures  3  and  4  illustrate  well 
the  stabilizing  influence  of  current  smoothing  and  oblate¬ 
ness.  However,  we  should  emphasize  that  the  radial  wall 
was  touching  the  separatrix,  that  the  separatrix  had  only 
modest  oblateness,  and  that  line  tying  was  present,  i.e., 
the  exterior  region  was  plasma.  In  the  next  section  we 
present  studies  of  stability  of  the  tilt  mode  with  walls 
further  removed  with  and  without  line  tying. 

When  our  equilibria  are  made  more  oblate  (i.e., 
larger  |s|)  than  those  shown  in  Figure  4,  a  new  instability, 
the  radial  shift  mode,  appears. 15  From  Figure  5,  we  see 
that  the  real  part  (i.e.,  cos  S  Fourier  component)  of  the 
poloidal  disolacement  $  is  roughly  a  rigid  radial  shift 

of  the  plasma.  The  imaginary  part  of  C  (the  sin  6  com- 

*\.P 

ponent)  is  again  tilt-like.  One  way  to  understand  the 
appearance  of  this  instability  for  more  oblate  plasmas  is 
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Fig.  5  —  Poloidal  displacement  ^  for  a  shift  eigenmode.  The  real  (cot  8 ) 
component  shown  in  (a)  shows  again  the  effect  of  magnetic  pressure  imbal¬ 
ance  near  the  perturbed  X-points.  The  imaginary  (sin  8 )  component  shown 
in  (b)  is  somewhat  tilt-like. 
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to  note  that  oblate  plasmas  are  produced  from  boundary 

conditions  (6)  or  (7)  by  increasing  the  amount  of  (negative) 

flux  passing  through  the  radial  wall,  i.e.,  by  increasing 
2  2 

5  s/oZ  at  r  =  a.  This  is  achieved  by  decreasing  e  in  (6) 

2  2 

or  by ‘  increasing  c^  in  (7).  Thus,  for  3  >/3Z  >  0,  the 

Grad-Shaf ranov  equation  (3b)  with  g('>)  =  0  implies 

r  3/3r  £{l/r)  3i(i/3r]<  0  or  3Bz/3r  <  0,  that  is  that  the  field 

index  is  negative.  On  this  basis,  we  might  expect  a  radial 

instability  like  the  radial  precession  mode  in  Astron  type 

geometries.^-®  The  simplest  model  for  this  instability 

applies  to  a  rigid  ring  of  current  so  weak  that  its  self 

forces  are  negligible.  For  this  model  instability  occurs 
2 

if  B  (Ne)  /4ttMc  I  +  v  <  0,  where  Ne  is  the  total  charge  of 
z 

the  current  carrying  species,  M  is  the  total  ring  mass,  I 
is  the  current  carried  by  the  ring  and  v  =  rB'/B  is  the 
(external)  field  index.  (The  first  term  in  this  inequality 
is  in  essence  a  finite  Larmor  radius  effect  and  may  be 
neglected  if  the  Larmor  radius  is  small  or  if  dissipation 
is  present.  ) 

A  related  approach  to  understanding  the  shift  insta¬ 
bility  is  to  note  that  very  oblate  equilibria  are  confined 
radially  by  strong  magnetic  fields  away  from  the  midplane 
z  =  0  but  by  rather  weak  fields  near  the  midplane,  so  that 
even  a  stable  rigid  shift  motion  may  have  a  weak  restoring 
force,  i.e.,  be  nearly  marginally  stable.  However,  such  a 
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shift  produces  magnetic  pressure  imbalance  at  the  perturbed 

X-points,  which  can  be  strong  enough  to  destabilize  the  mode. 

This  mechanism  is  closely  related  to  the  mechanism  by  which 

2 

the  external  tilt  in  oblate  plasmas  is  destabilized.  The 
inward  displacement  in  response  to  this  pressure  imbalan-  a 
is  clearly  evident  in  Figure  5. 


V. 


OPTIMUM  PARAMETERS  FOR  STABILITY 


In  this  section  we  present  results  of  a  systematic 
search  for  optimum  parameters  for  stability  to  both  the 
tilt  and  the  shift  mode.  We  have  performed  these  studies 
using  both  models  for  the  exterior  region,  namely  plasma 
exterior  (i.e.,  with  line  tying)  and  vacuum  exterior. 

In  Section  IV  we  studied  the  behavior  as  a  function 
of  5,  the  current  smoothing  parameter,  and  found  that 
stability  as  measured  by  the  position  of  the  axial  wall 
required  for  marginal  stability,  improved  monotonically 
with  6.  For  ^  =  0,  this  behavior  has  a  plateau  at  about 
5  =  .25. 

In  Figure  6  we  show  the  axial  wall  position  required 
for  marginal  stability  as  a  function  of  separatrix  shape  with 
for  o  =  .15  and  6  =  .25.  On  the  vertical  axis  is  the  ratio 
z^Zg,  where  zw  is  the  position  of  the  axial  wall  and  zg 
is  the  position  of  the  X-point,  where  the  separatrix  inter¬ 
sects  the  z-axis.  On  the  horizontal  axis  is  the  shape 
parameter,  the  ratio  of  half  length  to  radius  zQ/rs  (=1 
for  a  spherical  separatrix) .  Both  sets  of  runs  were  done 
with  line  tying  and  with  the  position  of  the  radial  wall 
relative  to  the  radial  position  of  the  separatrix  rw/rs 
fixed  at  1.25.  Note  that  in  both  cases  stability  is 
optimal  for  zs/rs%0.6.  The  equilibrium  with  5  =  .25  is 
clearly  more  stable  than  that  with  i  =  .15.  We  have  also 
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done  a  set  of  runs  with  <5  =  .35,  and  have  found  that  the 
marginal  stability  points  with  line  tying  lie  very  close 
to  the  6  =  .25  plot  shown  in  Figure  6. 

Figure  7  shows  the  stability  for  ^  =  0,  and  <5  =  .25 
both  with  and  without  line  tying.  The  plot  with  line  tying 
is  just  the  same  as  the  6  =  . 25  curve  shown  in  Figure  6. 
Without  line  tying,  stability  is  again  optimal  at  roughly 
zs/rs  ~  0.6.  Clearly,  line  tying  greatly  improves  the 
stability.  In  this  case,  with  our  fixed  radial  wall,  the 
axial  wall  is  permitted  to  be  about  nine  times  further  from 
the  separatrix  with  line  tying  than  without. 

For  Figure  8  we  have  taken  the  equilibrium  corre¬ 
sponding  to  optimal  stability  with  line  tying  in  Figure  7 
(rQ  -  .99999  and  e  =  -8.2  in  equation  6),  and  have  deter¬ 
mined  the  radial  wall  position  required  for  marginal 
stability  as  a  function  of  the  axial  wall  position.  Note 
again  that  line  tying  greatly  improves  the  stability. 

These  results  indicate  that,  in  the  presence  of  line  tying, 
the  radial  wall  alone  may  be  sufficient  to  stabilize  the 
n  =  1  modes.  Our  code  does  not,  of  course,  allow  us  to 
verify  this  directly,  because  we  require  a  finite  number 
of  axial  grid  points.  In  any  case,  it  is  clear  that  line 
tying  does  allow  us  to  remove  the  axial  wall  to  a  large 
distance,  as  required  by  the  moving-ring  reactor  scheme.^ 
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Fig.  7  —  Axial  position  of  the  wall  required  for  stability  as  a  function  of 
elongation  with  and  without  line  tying;  ^  ■  0,  e  -  0.25,  and  rw/r.  *  1.25. 


The  results  shown  in  Figure  8  show  clearly  that  when  the 
axial  wall  is  removed  to  a  large  distance,  line  tying  in 
the  radial  wall  is  responsible  for  stability. 

Finally,  we  have  looked  at  the  effect  on  stability 
of  a  flux  hole.  Our  results  with  current  smoothing  near 
the  separatrix  suggest  that  the  introduction  of  a  flux 
hole  should  be  stabilizing,  and  we  find  that  that  is  indeed 
the  case.  Again,  our  results  make  sense  if  we  recall  that 
the  displacement  near  the  X-point  plays  an  important  role 
in  destabilizing  both  the  tilt  and  the  shift  mode. 

Figure  9  shows  axial  wall  position  required  for 
marginal  stability  plotted  against  separatrix  shape,  with 
^h  =  ^ '  S  =  .05,  and  r^r  -  1.4.  The  optimally  shaped 

equilibrium  is  somewhat  more  stable  than  in  Figure  8,  where 
+  5  has  the  same  value  but  ^  =  0.  Now,  because  the  sur¬ 
face  of  our  current  carrying  equilibrium  plasma  lies  some¬ 
what  inside  the  separatrix,  the  wall  is  actually-  somewhat 
further  away  from  the  plasma  than  indicated  by  our  plot. 

This  can  be  seen  in  Figure  1,  which  shows  the  equilibrium 
corresponding  to  optimal  stability  with  line  tying  in 
Figure  9  (C-^  =  0.25  and  C ^  =  6.0  in  equation  7)  .  Equilib¬ 
rium  currents  flow  only  in  the  shaded  region,  with  the 
dashed  lines  showing  the  position  of  the  walls  at  marginal 
stability.  Once  again,  line  tying  greatly  improves  the 
stability.  The  results  of  our  runs  for  these  equilibria 
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Fig.  8  —  Axial  position  of  wall  vs.  radial  position  at  marginal  stability,  both 
scaled  to  the  relevant  separatrix  dimension;  6  =  0.25,  e  *  -8.2  (z,/rg  *  0.59). 
Results  are  shown  with  and  without  line  tying  of  the  open  field  lines. 


without  line  tying  are  consistent  with  the  results  obtained 
by  others  for  similar  equilibria . ^ 0 ' 11 ' ^  (For  comparison 
with  the  results  of  other  authors,  we  caution  that  the  wall 
radius  is  sometimes  expressed  as  the  distance  of  the  wall 
to  the  plasma,  i1  =  relative  to  the  minor  plasma  radius. 

In  this  normalization,  the  radial  wall  of  Figure  9  is  a 
distance  of  approximately  1.2  from  the  plasma,  and  the 
maximum  distance  of  the  axial  wall  without  line  tying  is 
about  .8.) 

Note  that  in  Figure  8,  with  r  /r  =  1.4,  the  marginal 

w  s 

position  of  the  axial  wall  for  the  equilibrium  without  a 
flux  hole  is  given  by  z^Zg  £  1.2.  Comparing  with  Figure  9, 
we  see  that  the  flux  hole  allows  us  to  move  the  walls 
further  from  the  separatrix,  this  despite  the  fact  that 
the  surface  of  the  current  distribution  has  moved  inside 
the  separatrix. 


VI. 


SUMMARY 


Our  tool  for  studying  the  st  ibility  of  force  free 
spheromak  equilibria  to  the  current  driven  tilt  and  radial 
shift  modes  is  the  time  dependent  magnetohydrodynamic  code 
FFMHD .  This  code  integrates  the  perturbed  vector  potential 

6A(  =  E,  x  b)  in  time  rather  than  the  displacement  For 

-v  a,  -v  ■v 

this  reason  it  is  well  suited  to  the  treatment  tilt  and 
radial  shift  modes,  whose  displacement  can  be  singular  at  the 
magnetic  neutral  points  (X-points)  on  the  axis  of  symmetry. 
The  external  region  may  either  be  treated  as  a  vacuum  or 
as  a  conducting  plasma  with  zero  equilibrium  current. 
Equilibria  treated  by  the  latter  model  are  generally  more 
stable  due  to  line  tying  of  the  open  field  lines. 

We  find  that  nearly  spherical  "Taylor  equilibria," 
having  V  x  B  =  yB  with  u  uniform  require  a  tightly  fitting 

'Xj  % 

cylindrical  wall  for  stability,  even  with  line  tying.  For 
equilibria  whose  current  goes  to  zero  more  smoothly  near  the 
separatrix,  and  with  a  conducting  plasma  in  the  exterior 
region,  the  walls  may  be  removed  a  considerable  distance. 

The  optimum  plasma  shape  for  stability  to  both  the 
tilt  and  the  radial  shift  modes  has  ratio  of  half  length 
to  radius  0.6,  nearly  independent  of  other  equilibrium 
parameters.  With  line  tying  on  the  open  field  lines,  the 
radial  and  axial  walls  can  be  removed  from  the  plasma  by  at 
least  50  percent,  relative  to  the  total  radius  and  length 
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of  the  separatrix.  The  axial  wall  may  be  removed  to  at 
least  several  times  the  plasma  length  with  the  radial  wall 
at  a  distance  of  over  1.2  times  the  plasma  radius.  These 
results  are  summarized  in  Figure  8.  With  vacuum  on  the 
open  field  lines,  the  optimum  shape  is  about  the  same,  but 
the  walls  must  be  considerably  closer. 
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4. 


APPENDIX  A:  NUMERICAL  CODE 

The  numerical  code  FFMHD  integrates  (9)  throughout 

the  region  enclosed  by  the  walls  or  (9),  (15),  depending 

upon  whether  the  exterior  is  modeled  by  a  conducting 

plasma  or  a  vacuum  (possibly  with  a  small  matching  region) . 

The  finite  difference  equations  used  are  based  upon  the 

simplest  central  differences  for  derivatives  with  respect 

2  2  2 

to  t,  r,  and  z.  Errors  are  of  order  At  ,  Ar  ,  Az  .  The 

Courant  condition  which  one  might  expect  for  such  a  scheme, 

2  2 

when  B  /o  3  vA  =  const.,  is  At  <  min  (Ar,  iz)/vA>  Empirically, 
we  find  that  a  slightly  smaller  time  step  is  required  for 
numberical  stability.  For  the  alternate  normalization 
p  =  const.,  the  required  time  step  is  about  20  percent 
smaller  (where  vA  now  represents  the  maximum  Alfven  speed) . 

As  discussed  in  Section  III,  we  expect  that  a  model  with 
low  density  in  the  exterior  region  [in  order  to  obtain 
correct  growth  rates  with  (16)]  would  require  extremely 
short  time  steps. 

The  boundary  conditions  on  the  perturbed  vector 
potential  6 A  at  the  conducting  walls  have  been  discussed 
in  Section  III.  Boundary  conditions  are  also  required  on 
the  symmetry  axis  r  =  0  because  of  the  coordinate 
singularity. 
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We  describe  the  boundary  conditions  at  r  *  0  for 
each  toroidal  mode  number  n  separately.  Assuming  the 

components  5Ax#  SA^  to  be  analytic  in  x,  y,  z,  i.e., 

expanding  for  fixed  z 

5  A  =  A  +  Bx  +  Cy  +  •  ■  • , 

x 

£Ay  =  D  +  Ex  +  Fy  +  •••,  (Al) 

6Az  =  H  +  lx  +  Jy  +  •  •  •  , 

we  find 

6Af  =  ^(B  +  F) r  +  A  cos  6  +  D  sin  6 

+  ^(B  -  F) r  cos  26  +  |(C  +  E) r  sin  26  +  •••, 

5Aq  =  ^(E  -  C)r  +  D  cos  9  -  A  sin  6  (A2) 

+  ^-(E  +  C)r  cos  29  +  ^(F-B)r.sin  29  +  •••, 


oA  =  H  +  I  r  cos  6  +  J  r  sin  9  +  • • • . 
z 

As  in  Section  III,  we  write  £A  =  Re[(5Ar  +  icA1)exp  (in9)l 

•x,  x  x  c 

*  5Ar  cos  n  9  -  5A1  sin  n  9.  For  n  =  0  we  conclude 
%  % 

r*  j 

A  =  0  =  0  and  therefore  i A  =  £Aa  =  0.  (Imaginary  parts 

it  y 


are  automatically  zero.)  Also,  since  the  neglected  terms 

of  6A  are  O(r^),  we  have  (3/9r)5A^  =  0. 
z  z 


For  n  *  1,  we  find  B  +  F=*B-F=G,  C  +  E»C-E« 

0,  so  that  (3/5r)5A  *  (j/dr)6Aa  *  0  for  real  and  imaginary 

r  u 

parts  and  6Ar  *  5A~ ,  ?a 1  *  -6A^ .  We  also  find  6A*  ■  A1  *  0. 

F  r  u  r  6  z  z 

2 

For  n  ■  2  we  find  5A  ^  dA.  v  r  as  r  •*  0  and  6A  a,  r 

r  9  z 

as  r  -*•  0.  It  is  easy  to  see  that  for  n  _>  2  this  generalizes 

to  6A  %  6Afl  >.  r11-1,  6A  ^  rn. 
r  y  z 

The  conditions  for  n  *  0  are  easily  implemented  in 

the  code.  Similarly,  for  n  ^  2,  it  is  adequate  to  set  all 

three  components  of  dA  equal  to  zero  at  r  =  0.  For  n  =  1 

the  issue  is  complicated  by  the  fact  that  the  conditions 

given  overdetermine  <5Ar  and  cA^  at  r  *  0,  and  implementing 

any  subset  of  these  conditions  leads  to  numerical  instability. 

We  were  able  to  overcome  these  difficulties  by  setting  dA^ 

and  5Ag  both  equal  to  the  average  of  the  values  of  5A^  and 

c'aJ  obtained  by  using  the  conditions  (5/ dr)  5Ar  *  (3/jr)  dA^  *  0. 
y  r  y 

X  -  IT 

The  analogous  procedure  was  used  for  5Ar  and  iA^. 

When  the  code  is  run  with  the  option  that  initial 
conditions  are  generated  in  the  code,  we  insure  that  the 
initial  conditions  obey  the  boundary  conditions  at  r  =»  0 
as  well  as  those  at  the  conducting  wall.  (The  other  option, 
which  speeds  up  convergence  to  the  most  unstable  mode,  is 
to  input  the  vector  potential  from  a  previous  run.) 

The  typical  grid  we  use  to  integrate  the  finite 
difference  form  of  (9)  or  (9),  (15)  has  40  or  50  points 
radially  and  100  points  axially.  We  have  performed 


I 


convergence  tests  on  finer  grids  (up  to  70  x  170)  and  con- 

% 

elude  that  the  relative  grid  error  in  y  *  (-dW/T)  is 

n 

less  than  three  percent.  In  the  same  set  of  runs  the 
improvement  in  the  accuract  of  y  (see  Section  III)  was 
much  more  pronounced.  In  the  transient  phase,  i.e.,  for 
short  times,  energy  is  conserved  in  the  code  to  within  a 
few  percent  (the  change  in  total  energy  relative  to  the 
change  in  kinetic  energy)  for  a  40  x  100  grid  and  to  with¬ 
in  a  fraction  of  a  percent  for  a  70  x  170  grid.  By  the 
time  an  unstable  normal  mode  has  evolved,  energy  conser¬ 
vation  is  worse.  We  trace  this  behavior  to  singular 
behavior  near  the  X-points.  However,  in  no  cases  that  we 
have  tested  does  this  reflect  itself  in  poor  accuracy  for 
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APPENDIX  B:  EXTERIOR  REGION  MODELS 


Our  purpose  here  is  to  explain  in  more  detail  various 
properties  of  our  modified  magnetohydrodynamic  equations 
in  the  exterior  region.  Specifically,  we  prove  that  normal 
modes  of  (9), (15)  have  continuous  tangential  components  of 
<5 A,  and  that  the  model  with  plasma  exterior  is  more  stable 
than  the  one  with  vacuum  exterior  because  of  line  tying. 

We  first  prove  that  regular  (i.e.,  bounded)  normal 
modes  of  the  system  with  a  vacuum  exterior  described  by 
(9),  (15)  have  continuous  tangential  components  n  x  SA  at 
the  separatrix,  even  when  the  exponential  factor  in  (15) 
is  zero.  We  identify  three  components  of  the  perturbed 
vector  potential  at  the  separatrix  £A„,  SA,. ,  and  SA. _. 
These  are,  respectively,  the  normal  component,  the  tangen¬ 
tial  component  parallel  to  B,  and  the  tangential  component 

perpendicular  to  B.  We  also  split  the  curl  operator  into 

% 

tangential  and  normal  components,  so  that  SB  =  7  x  SA  = 

%  'V 

x  {A  + 

:  -vn 


*®n  + 

4St' 

where  SB  =  7^  x  5 A. ,  and  5B.  = 
^n  t  --t  %t 

7n  x 

c£t* 

We 

also  find  5j  =  7  x  SB  =  5j 

<\T  'v  <n 

'Jn 

7fc  x 

7t 

x  SA  +  7.  x  7  x  SA  and  Sj 

<\,n  t  n  <\,t  \ 

7n  x 

7t  X 

SA 

'vn 

+  7  x  7  x  SA. .  Therefore 
n  n  %t 

equations  (9) , (15) ,  with  B  /p  constant  for  convenience, 
and  near  the  separatrix  where  u(^)  can  be  neglected,  are 
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Y  2<5  A. 


y2s^t2 


5A.  ,  =  0  for  ip  <  0 
xtl 


Y  ^  '  A 

Y  tl 


•  2*tl 


for  p  >  0 


Now,  for  contradiction,  assume  5A,-  is  discontinuous  at 

'vtz 

i>  -  0.  Then,  from  (B2)  ,  either  -5A.  ~  must  behave  like 
5'(p)  or  6An  must  behave  as  5  {•;•).  Both  of  these  conclusions 
contradict  the  assumption  that  5 A  is  bounded,  so  that  6A. _ 
is  continuous.  Next,  assume  that  <SA.,  is  discontinuous  at 
p  *  0.  Then,  from  (Bl)  and  since  5Afc2  is  continuous,  we 
see  that  5An  must  behave  as  5(d),  which  is  again  a  contra- 

v 

diction.  Therefore  the  tangential  components  of  5A  are 
continuous  at  p  =  0.  Incidentally,  a  similar  argument 
cannot  be  applied  to  prove  the  continuity  of  the  normal 
component  of  cA,  because  the  only  normal  derivative  of 
5A„,  in  5j^,  is  balanced  by  a  second  normal  derivative  of 

'\,t  - 

5Afc.  The  fact  that  the  system  (9),  (15)  requires  continuity 
of  5Afc  but  allows  discontinuity  of  5An  explains  why  FrMHD 
has  no  numerical  problems  when  run  in  the  exterior  vacuum 
mode,  even  when  the  exponential  factor  in  (15)  is  so  small 
that  it  is  essentially  zero  on  every  grid  point  in  the 
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exterior  region.  Also,  it  is  easy  to  see  that  the  above 

2 

proof  does  not  require  that  B  /o  be  constant,  or  even  that 

p  be  continuous  at  'i>  =  0.  (Note  that  if  we  pick  the  fic- 

2  2 

titious  density  p  so  that  vA  =  c  in  the  exterior  region, 
then  (15)  merely  states  that  7  x  6B  is  balanced  by  the 
displacement  current  in  the  exterior,  and  therefore  that 
the  perturbed  current  is  zero.) 

In  the  previous  paragraph  we  have  completed  the 
proof  that  the  system  (9) , (15)  has  identical  magnetohydro¬ 
dynamic  stability  properties  as  the  physical  system  with 
vacuum  on  the  open  field  lines,  since  it  has  the  same  con¬ 
tributions  to  5W  in  plasma  and  in  vacuum  [c.f.  (11)]  and 

since  the  usual  boundary  condition  ft  x  5A  =  n  x  (£xB)  = 

-  ft  •  5  B  holds.  We  now  turn  to  a  proof  that  the  system 
%  % 

described  by  alternate  equations  of  motion  (9) ,  with  a 
plasma  on  the  open  field  lines,  is  more  stable  than  the 
system  with  vacuum  exterior.  This  stabilizing  influence 
is  due  to  the  fact  that  we  have  a  normal  component  of  the 
equilibrium  field  B  at  the  walls,  i.e.,  we  have  line  tying. 
Indeed,  it  is  shown  in  Ref.  8  that  if  ft  •  B  =  0  at  the 
walls,  and  if  no  mode  rational  surface  exists  in  the  plasma, 
the  condition  B  •  5A  =  0  can  always  be  satisfied  by  a  gauge 
transformation . 
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The  remaining  question  is  whether  the  additional  con¬ 
straints  on  the  vector  p itential  in  the  exterior  plasma 

model,  namely  B  •  6 A  =  0  and  continuity  of  6A  at  the 

%  % 

separatrix,  are  physical  constraints  (i.e.,  constraints 
on  the  fields)  or  whether  they  can  be  imposed  by  a  gauge 
transformation.  If  we  try  to  impose  the  former  condition, 

i.e.,  try  to  satisfy  B  •  7x  +  B  •  6A  =  0,  we  find 

%  %  v 


X 


dSL 

% 


(B4) 


where  the  integral  is  along  a  field  line  and  C  is  the  value 
of  x  at  the  wall  (this  must  be  a  constant  by  the  condition 
that  the  tangential  component  of  <$A  in  both  gauges  must 
vanish  at  the  wall) .  Therefore,  for  any  field  line  that 
passes  through  the  walls  twice  (see  Figure  1)  we  must  have 


/*  • 


d£  =  0.  For  n  =  0  (axisymmetric)  modes  this  means 


/cB 

* 


that  the  perturbed  toroidal  flux  $  cB  •  ndA  between  any 
two  flux  surfaces  p  =  const,  in  the  exterior  region  must 
be  zero.  For  modes  with  n  /  0  we  have,  from  Appendix  A, 
£A^  =  0  on  the  axis  r  =  0.  Therefore  the  perturbed 
toroidal  flux  enclosed  by  the  z-axis  and  any  flux  surface 
p  -  const,  that  intersects  the  walls  must  be  zero.  There¬ 
fore  the  constraint  B  •  5A  =  0  is  a  physical  constraint 
and  has  a  stabilizing  influence.  As  an  aside,  we  comment 


that  any  closed  flux  surfaces  in  the  flux  hole  p  >  p. 
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[c.f.  equation  (4)]  may  be  considered  to  be  in  vacuum  or 
to  contain  plasma.  For  such  flux  surfaces,  with  Bg  =  0, 
the  condition  B  •  cA  =  0  can  be  imposed  by  a  gauge  trans- 
formation  only  if  the  integral  over  the  closed  field  line 


/ 


5 A  •  dll  =  0.  Since  this  implies  that  the  perturbed 
x  % 


toroidal  flux  is  zero  within  the  flux  surface,  having 
plasma  on  the  closed  field  lines  in  the  flux  hole  is  also 
stabilizing.  (This  is  a  special  case  of  the  theorem  of 
Reference  8,  that  the  constraint  ^  =  0  is  stabilizing 

if  a  mode  rational  surface — here  q  =  0 — exists.) 

Finally,  the  boundary  condition  on  6^  is  a  further 
physical  constraint.  Indeed,  let  us  require  that  the  same 
gauge  function  of  (B.4)  satisfy  3x/^n  =  ( 5Ap  -  6A •  n 
at  =  0,  where  3x/3n  is  the  normal  derivative,  6a  =  £  x  b 

'Vp  X  x 

is  the  vector  potential  inside  the  plasma  ( <  0)  and  oA 


is  the  vacuum  vector  potential.  We  find  from  (B4) 


fe 


<5A  •  di  =  -  An  •  6A  , 

a.  'v,  %  P 


'.V 


( B5 ) 


=  A  I  B  I  £ 
a  1 


The  integral  on  the  left  is  over  a  curve  consisting  of  the 
z-axis  exterior  to  the  plasma  and  part  of  the  separatrix 
up  to  the  point  in  question,  and  returning  on  a  nearby 
field  line,  separated  by  a  distance  A  at  the  terminal  point 
on  the  separatrix.  Thus,  the  continuity  boundary  condition 
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on  5An  imposes  a  physical  constraint  on  the  flux  through 
the  area  bounded  by  this  curve,  and  is  therefore  also 
stabilizing. 
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